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Abstract 

This paper provides a new algorithm for solving inverse problems, based on the 
■ minimization of the norm and on the control of the Total Variation. It consists in 

, relaxing the role of the Total Variation in the classical Total Variation minimization 

I approach, which permits us to get better approximation to the inverse problems. 

The numerical results on the deconvolution problem show that our method outper- 
forms some previous ones. 
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O ■ 1 Introduction 

Environmental effects and imperfections of image acquisition devices tend to degrade the 
quality of imagery data, thereby making the problem of image restoration an important 
part of modern imaging sciences. Very often the images are distorted by some linear 
transformations. In this paper, we consider the reconstruction of the unknown function 
5^ \ f from the following inverse problem: 

Y{x) = Kf{x) + e(x), x G / = {0, iV - 1}\ (1) 

where Y is the observed data, /C is a continuous linear operator from R^^ to i?*^, / is the 
original image, e is a white Gaussian noise of mean and variance a^, > 1 and M > 1 
are integers. We want to recover the original image / starting from the observed one Y . 
When the operator /C is formulated as a convolution with a kernel, this inverse problem 
reduces to the image deconvolution model in the presence of noise. Of particular interest 
is the case where the operator K is the convolution with the Gaussian kernel: 

K:f{x) = g*f{x) = Y,9{^-y)f{y)^ ^e/, (2) 
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where 



and (T;, > is the standard deviation parameter. This type of inverse problems appears 
in medical imaging, astronomical and laser imaging, microscopy, remote sensing, photog- 
raphy, etc. 

The problem ([1]) can be decomposed into two sub-problems: the denoising problem 

Y[x) = u[x)^e[x), x^I (4) 

and the possibly ill-posed inverse problem 

u [x] = )Cf (x) , X E I. (5) 

Thus, the restoration of the signal / from the model ([T]) will be performed into two steps: 
firstly, remove the Gaussian noise in the model (jl]); secondly, find a solution to the problem 
©. 

For the first step, there are many efficient methods, see e.g. Buades, Coll and Morel 
(2005 [2]), Kervrann (2006 [3), Lou, Zhang, Osher and Bertozzi (2010 Q), Polzehl 
and Spokoiny (2006 jl7|), Garnett, Huegerich and Chui (2005 0), Cai, Chan, Nikolova 
(2008 fl), Katkovnik, Foi, Egiazarian, and Astola ( 2010jl3l), Dabov, Foi, Katkovnik 
and Egiazarian (2006 j^), and Jin, Grama and Liu(2011 

The second step consists in finding a convenient function / satisfying (exactly or 
approximately) 

JCf-u = 0, (6) 

where u is obtained by the denoising algorithm of the first step. To find a solution to 
the ill-posed problem ([HD, one usually considers some constraints which are believed to be 
satisfied by natural images. Such constraints may include, in particular, the minimization 
of the norm of / or V/ and the minimization of the Total Variation. 

The famous Tikhonov regularization method (1977 [19|] and 1996 j^) consists in solving 
the following minimization problem: 

f{x) = argmin ||7||^ (7) 

subject to ^ 

\\]Cf-u\\l = 0. (8) 

Using the usual Lagrange multipliers approach, this leads to 

/ = argmin ||7||2 + 1a||/c7- «||^, 
/ ^ 

where A > is interpreted as a regularization parameter. 

An alternative regularization (1977 [lt|) is to replace the norm in ([7]) with the 
semi-norm, which leads to the well-known Wiener filter for image deconvolution: 

/ = argmin ||V7||2 (9) 
/ 
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subject to ([8]), or, again by the Lagrange multipliers approach, 



/ = argniin || V/||^ + ^A||/C/ - (10) 

where A > is a regularization parameter. 

The Total Variation (TV) regularization called ROF model (see Rudin, Osher and 



Fatemi [l8|), is a modification of with ||V/||2 replaced by ||V/||] 



/ = argmin ||V/||i (11) 
/ 

subject to ([8]), or, again by the Lagrange multipliers approach, 

/ = argmin l|v7||i + ^A||/c7- u\\l. (12) 

where A > is a regularization parameter. 

Due to its virtue of preserving edges, it is widely used in image processing, such as 



blind deconvolution p, [lOl, ll5| , inpaintings j5(| and super- resolution [161 . Luo et al [14 
composed TV-based methods with Non-Local means approach to preserve fine structures, 
details and textures. 

The classical approach to resolve f|T2|) . for a fixed parameter < A < +oo, consists in 
searching for a critical point / characterized by 

L(7) + A/C*(/c7-i/) = 0, (13) 

where L{f) = V^jj, |V/| being the Euclidean norm of the gradient V/. The usual 

technique to find a solution / of f[T^ is the gradient descent approach, which employs the 
iteration procedure 

fk+i = fk-h (Lfk + A/C*(/C/fc - u)) , (14) 

where /i > is a fixed parameter. If the algorithm converges, say fk — /, then / satisfies 
the equation (1131) . However, this algorithm cannot resolve exactly the initial deconvolution 
problem Qj, for, usually, L{f) ^ 0, so that, by f|T3|) . /C/ — m 7^ 0. This will be shown in 
the next sections by simulation results. 

In this paper, we propose an improvement of the algorithm f[T^ . which permits us to 
obtain an exact solution of the deconvolution problem ([6]) when the solution exists, and 
the closest approximation to the original image when the exact solution to ([6]) does not 
exist. ^ 

Precisely, we shall search for a solution / of the equation 

]C*{]Cf-u) = 0, 

with a small enough TV. Notice that the last equation characterizes the critical points 
/ of the minimization problem minj ||/C/ — 'u||2 which, in turn, whose minimiser gives 
the closest approximation to the deconvolution problem iQ, even if it does not have any 
solution. In the proposed algorithm, we do not search for the minimization of the TV. 
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Instead of this, we introduce a factor in the iteration process (IT^ to have a control on 
the evolution of the TV. Compared with the usual TV minimization approach, the new 
algorithm retains the image details with higher fidelity, as the TV of the restored image 
is closer to the TV of the original image. 

Our experimental results confirm this point, and show that our new algorithm out- 
performs the recently ones proposed in the literature that are based on the minimization 
problems ([ID]) and ([12]). 

2 A new algorithm for the inverse problem 

2.1 Controlled Total Variation regular izat ion algorithm 

We would hke to find / satisfying the deconvolution equation ([6]) with a small Total 
Variation. Since an exact solution to ([6]) may not exist, instead of the initial problem we 
consider the minimization problem 

min \\ICf-u\\l (15) 
/ 

This leads us to search for critical points / characterized by the equation 

]C*{]Cf-u)=0, (16) 

whose TV is small enough. As explained in the introduction, the usual algorithm f lT4|) 
(classic TV) does not lead to a solution of f lT6]) if it exists. In fact, we have observed 
by simulation results that usually ||L/fc||i does not converge to with k — )■ +oo, so that 
lC*{K,fk — u) does not tend to neither. 

For example taking "Lena" as test image blurred by the Gaussian convolution kernel 
([3]) with standard deviation cr;, = 2, the value of increases (see Figured) (b), where 

the classical TV curve correspond to the dashed line) as log/c increases. On the other 
hand. Figure [T](c) shows that the value of ||/C*(/C/fc — m)||i decreases when log/c increases 
in the interval [0,3], but it remains bounded from below by — 1 ~ 2980 and does not 
decrease towards when log/c varies in the interval [3,8]. 

These considerations lead us to the following modified version of f[T4|) : 

fk+i = fk~h {OkLfk + A/C*(/C/,. - u)) , (17) 

where h > 0, A > are as in f[T^ and 0^ > is a sequence decreasing to 0. We shall use 
6k = 9'', for some 9 G (0, 1). This ensures that, if the algorithm converges, say fk — )■ /, 
then 

}C*{}Cf-u)=0. (18) 

Notice that, compared with the initial algorithm ([T^ . we put a small weight 9k to 
the TV in the minimization problem ([12]) . thus diminishing the role of this term in the 
minimization, in order to get an exact solution to the ill-posed problem fll8p . Here we do 
not search for the minimization of the TV, but just for a control on the TV. The control 
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of TV rather than the minimization of TV permets us to better preserve image edges and 
details. 

The new algorithm given by fll7p will be called Controlled Total Variation regular- 
ization algorithm (CTV). 

We show by simulations that the new algorithm improves the quality of image restau- 
ration. In Figures [1] and O we see that the TV of the restored image by the proposed 
CTV algorithm is larger than that of the restored image by the classical TV minimiza- 
tion approach, but closer to the TV of the original image. This also shows that in the 
classical TV minimization algorithm the TV of the restored image is too small and that 
the classical approach may not give the best solution. 

On the other hand, by simulation we have also seen that the TV in the iterative process 
( IT7|) cannot be suppressed completely, i.e. we cannot take O^. = 0. 

In the case of the Gaussian blur kernel. Figured] (c) shows that values of \ \}C*{}Cfk—u)\\i 
computed by the new algorithm CTV and those by the classical TV algorithm are nearly 
the same, when log k varies in the interval [0, 3]. However, for log k varying in the interval 
[3, 8], the values of ||/C*(/C/fc — u)||i by the classical TV algorithm remain constant, while 
those computed by our algorithm continue to decrease to 0. Accordingly, we see in Figure 
[T] (d) that the PSNR (the definition of PSNR see Section 13.11) values by the classical TV 
algorithm and by the new proposed one, are almost the same for log k varying in the 
interval [0,3], imitating the evolution of ||/C*(/C/fc — however, for log A; varying in the 
interval [3, 8], the PSNR value by the classical TV algorithm does not increase any more, 
while the corresponding PSNR by the new algorithm continues to rise. In the case of 9 x 9 
box average blur kernel. Figure |2] displays similar evolutions as those of the Gaussian blur 
kernel case. Figure E] shows the images restored by our method and the classical TV 
method respectively. Our method keeps much more details of the image than the classical 
TV method. 

2.2 Direct Gradient Descent 

Suppose that K, a linear operator from R'^^ to R^^ . To solve the deconvolution problem 
dnj, inspired by ( 1T7|) . we can also consider the following modified version, that we call 
Direct Gradient Descent (DGD): 

fk+i = fk-h {OkLfk + A(/C/fc - u)) , (19) 

where < — 0. If the algorithm converges, say fk ^ f, then / satisfies exactly the 
equation (|6]). Therefore we could expect that this algorithm gives a better solution to 
the deconvolution problem. By simulations, this is shown to be the case when K, is the 
Gaussian blur kernel and the noise is absent. In this case we find that the DGD approach 
restores the blurred image more efficiently than classical TV and the CTV regularization 
approaches introduced above. The results of comparison between the classical TV, CTV 
and DGD are displayed in Figures [1] and |H 

However, our simulations also show that the algorithm DGD is very sensible to the 
noise and often may not converge. For example, when /C is the Gaussian blur kernel, but 
in the presence of the noise, the method does not give good restoration results; when /C 
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Classical TV and CTV (Gaussian blur kernel) 
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Figure 1: Comparison between classical TV and CTV. The test image "Lena" is convoluted 
with Gaussian blur kernel of standard deviation o"b = 1. We choose h = 0.1, A = 1, max Iter = 
3000 and 9 = 0.98. (a) The evolution of TV{ fk) as a function of log A;, (b) The evolution of 
log(| |L/fc| |i + 1) as a function of log k. (c) The evolution of log(||/C*(/C/fc — ii)||i + 1) as a function 
of log k, (d) The evolution of PSNR value as a function of log k. 
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Classical TV and CTV (box average blue kernel) 
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Figure 2: Comparison between classical TV and CTV. The test image "Lena" is convoluted 
with 9x9 box average blur kernel. We choose h = 0.1, A = 1, maxlter = 3000 and 6 = 0.998. 
(a) The evolution of TV{fk) as a function of \ogk. (b) The evolution of log(||L/fc||i + 1) as a 
function of log A:, (c) The evolution of log(||/C*(/C/fc — n)||i + 1) as a function of log A;, (d) The 
evolution of PSNR value as a function of log k. 
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Restoration results with classical TV and CTV 




(c) CTV, PSNR= 32.06db" (d) classical TV, PSNR= 28.33db 



Figure 3: (a) and (b) are the images restored from the test image "Lena" convoluted with 
Gaussian blur kernel with (T5 = 1 by our method and by classical TV method respectively, (c) 
and (d) are the images restored from the test image "Lena" convoluted with 9x9 box average 
blur kernel by our method and by classical TV method respectively. 
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is the 9x9 box average kernel, the algorithm diverges and the restoration results are not 
good with or without noise (see Figure E]). 

Notice also that the CTV regularization approach applies even if the linear operator /C 
is from to i?*^, where M is possibly different from A^^, contrary to the DGD approach 
which applies only when M = N"^ . 

The aforementioned arguments show that the CTV regularization approach f|T71) should 
be proffered to the Direct Gradient Descent algorithm f[T^ . especially when the blurring 
is accompanied with a noise contamination. 



3 Numerical simulation 
3.1 Computational algorithm 

In this section we present some numerical results to compare different deblurring methods. 
We shall use CTV regularization algorithm to restore several standard blurred and noisy 
images and compare it to the Tikhonov + NL/H^ (Tikhonov regularization followed by 
a Non-Local weighted approach) and Tikhonov + NL/TV (Tikhonov regularization 
followed by a Non-Local weighted TV approach) algorithms presented in When no 
noise is present, we shall also give simulation results by the DGD approach. 

Algorithm : BM3D + CTV(OPF + CTV) 

Input: Data Y and the degradation model /C. 
Step 1. Denoising. 



to 



Remove Gaussian noise from Y by BM3D [7] or by Optimization Weight Filter [11 
get the denoised image u. 
Step 2. Deconvolution. 

Set e = 5e~^ and choose parameters A, h, 6, k and maxlter = 3000. 

Initialization: k = and /o = u. 

(a) compute fk+i by the gradient descent equation f|T9|) 

(b) if ||/C*(/C/fc+i - u) - }C*{)Cfk - u)\\2> e and k < maxlter set k = k + 1 and go 
to (a). 

Output: final solution /fe+i. 

Note that If there is no noise, we just do Step 2. 

In the algorithm we use the control condition ||/C*(/C/fc+i — u) — fC*{lCfk — u)\\2 < e 
since we would like to ensure that ]C*{]Cfk — -u) — )■ 0. We could replace this condition by 
\\fk+i-fkh<e. 

The performance of the CTV regularization algorithm is measured by the usual Peak 
Signal-to-Noise Ratio (PSNR) in decibels (db) defined as 

PSNR = 10 logio 



MSE' 

where 



with / the original image and / the estimated one. 
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Figure 4: Comparison between CTV and DGD. Tlie test image "Lena" is convoluted with 
Gaussian blur kernel of standard deviation cTf, = 1. We choose h = 0.1, A = 1, max Iter = 
3000 and 9 = 0.98. (a) The evolution of TV{fk) as a function of log A;, (b) The evolution of 
log(||L/fc||i + 1) as a function of log/c. (c) The evolution of — fk\\i as a function of log/c, 

(d) The evolution of PSNR value as a function of log k. 



Table 1: The statistics of blur and noise we add to the images 



Image 


Blur kernel 




Gaussian noise 1 


Gaussian noise 2 


Gaussian noise 3 


Shape 


Gaussian with ai, = 


2 


cr„ = 


(7n = 10 


an = 20 


Lena 


Gaussian with o-j, = 


1 


u„ = 


(T„ = 10 


an = 20 


Barbara 


Gaussian with ai, = 


1 


u„ = 


an = 5 


an = 20 


Carmeraman 


9x9 box average 




cr„ = 


an = 3 


an = 10 



10 



DGD with box average blur kenel 
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Figure 5: The test image "Lena" is convoluted with 9x9 box average blur kernel. We choose 
h = 0.1, A = 1, ma^ Iter = 3000 and d = 0.998. (a) The evolution of log{TV{fk) + 1) as a 
function of log/c. (b) The evolution of logdl/fc+i — + 1) as a function of log/c. 




(a) Shape 150 X 150 (b)Lena 256 X 256 (c) Barbara 256 X 256 (d) Cameraman 256 X 256 

Figure 6: The test images 
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Table 2: The PSNR values of the simulation results for different methods 



Image 


Blur kernel 




Noise 


Tikhonov 


Tikhonov 




riiVloJJ 












+ NL/Hi 


+ NL/TV 


+ CTV 


+ CTV 










Un = 


23.62db 


23.62db 


32.28db 


32.28db 


35.49db 


Shape 


Gaussian with ai, = 


2 


O-n = 10 


27.53db 


28.71db 


29.13db 


29.15db 


— .— db 








O-n = 20 


25.02db 


26.58db 


27.75db 


27.50db 


— .— db 








CT„ = 


37.80db 


37.80db 


42.05db 


42.05db 


53.53db 


Lena 


Gaussian with cri, 


1 


o-„ = 10 


27.56db 


27.55db 


29.08db 


29.53db 


— .— db 








O-n = 20 


25.29db 


25.26db 


27.48db 


27.84db 


— .— db 








CT„ = 


40.86db 


40.86db 


45.45db 


45.45db 


56.82db 


Barbara 


• with (T(, = 


1 


<Jn = 5 


28.17db 


28.25db 


29.48db 


30.35db 


— .— db 








O-n = 20 


23.31db 


23.36db 


24.67db 


25.45db 


— .— db 








CT„ = 


27.22db 


27.22db 


30.37db 


30.37db 


— .— db 


Carmeraman 


9x9 box average 




0"n = 3 


25.40db 


25.43db 


25.78db 


26.07db 


— .— db 








o-„ = 20 


22.10db 


22.28db 


22.58db 


22.80db 


— .— db 



We test all the methods on four images: a synthetic image, Lena, Barbara and Camera- 
man (see Figure [6]) with several kinds of blur kernels and several values of noise variances 
as listed in Tabled! The synthetic image is referred to as Shape since it contains geometric 
figures. The Cameraman image is of high contrast and has many edges, while the Lena 
and Barbara images contain more textures. 

3.2 Simulation results 

We use BM3D or Optimization Weights Filter to remove Gaussian noise. Then we use 
our CTV algorithm to inverse the denoised blurred image. For each method, we present 
the PSNR results along with the applied parameters. We list the PSNR values of the 
used methods in Table EJ The Direct Gradient Descent method inverse excellently the 
convoluted image with Gaussian blur kernel without any noise, and the PSNR values are 
much higher than PSNR values of the other algorithms. However, the method works well 
only for the Gaussian blur kernel without any noise. For the noisy and blurred images, 
the method BM3D + CTV is the best way to reconstruct these images. Figures [TJ El [9] 
and [To] display the reconstructed images restored by the different methods and the square 
of the difference between the restored and original images. Figures [TJ [8] and [10] show that 
the method BM3D + CTV is the best one to remove the noise and blur. From Figure 
in] we see that BM3D + CTV works very well in reconstructing the image blurred with 
Gaussian blur kernel without any noise, but the method DGD can work much better than 
BM3D + CTV. 

4 Conclusion 

The Total Variation (TV) approach is a classical method to reconstruct images from the 
noisy or blurred ones. We have reconsidered deeply this approach and propose a modified 
version which we call Controlled Total Variation (CTV) regularization. The proposed 
algorithm is based on a relaxation of the TV in the classical TV minimization approach 
which permits us to find the best approximation to the solution of inverse problems. 
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Original and degraded image Tikhonov + NL/H\ PSNR= 27.53db 




BM3D + CTV, PSNR= 29.irKll) 




Figure 7: A 150 x 150 image with Gaussian blur cr;, = 2 and Gaussian noise (7„ = 10, the 
reconstructed image by different methods and their square errors. 
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Original and degraded image 




Tikhonov + NL/TV, PSKR= 27.55clt) 




Figure 8: Lena 256 x 256 with Gaussian blur 
structed image by different metfiods and their 




OPF + CTV, PSNR= 29.08(11) 




CTfe = 1 and Gaussian noise (T^ — 10, the recon- 
square errors. 
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Figure 9: Barbara 256 x 256 with Gaussian blur aj, = 1 and Gaussian noise (7„ = 0, the 
reconstructed image by different methods and their square errors. 
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Figure 10: Caxmeraman 256 x 256 with box average kernel and Gaussian noise cr„ = 20, the 
reconstructed image by different methods and their square errors. 
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Numerical simulations show that our CTV algorithm gives superior restauration results 
compared to known methods. 
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